Surface melting of the vortex lattice 
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We discuss the effect of an (afe)-surface on the melting transition of the pancake-vortex lattice 
in a layered superconductor within a density functional theory approach. Both discontinuous and 
continuous surface melting are predicted for this system, although the latter scenario occupies the 
major part of the low-field phase diagram. The formation of a quasi-liquid layer below the bulk 
melting temperature inhibits the appearance of a superheated solid phase, yielding an asymmetric 
hysteretic behavior which has been seen in experiments. 

PACS numbers: 74.25.Qt, 74.25.Dw, 68.35.Rh 



The transition between vortex solid and vortex liq- 
uid phases observed in the mixed phase of high-Tc su- 
perconductors has generated renewed interest in the 
problem of melting. In common with other discon- 
tinuous phase transitions, melting should involve the 
appearance of both metastable undercooled liquid and 
overheated solid phases. Hence, hysteretic behavior is 
expected upon cycling through the transition. How- 
ever, experiments on the layered, high- Tc superconduc- 
tor BiSSCO reveal an asymmetric hysteresis, charac- 
terized by the appearance of only the supercooled liq- 
uid and no overheated solid. Similar behavior is dis- 
played by ordinary crystals 0, 0| , where such asymme- 
try is understood to be a consequence of surface (pre-) 
melting: surfaces act as nucleation centers for the liquid, 
thereby inhibiting the metastable solid above the melting 
transition. However, such surface melting is not generic 
and there are experimental systems where the surface re- 
mains solid up to the bulk melting transition 0|. In this 
letter, we study the effects of an (a6)-surface on vortex 
lattice melting, showing that as the strength of the mag- 
netic field is varied, the same surface may exhibit either 
'surface non-melting' or 'surface melting' behavior. The 
latter scenario applies to the major part of the low-field 
phase diagram, in agreement with experiments 0. 

Early studies |^] of simple crystals have focused on the 
solid phase and have demonstrated that the surface turns 
unstable before the bulk melts. Going beyond such a sta- 
bility analysis is a difficult task, as a theory is required 
which describes both solid and liquid phases in a unified 
manner. Qualitative insight can be gained from a Lan- 
dau theory 6] by including the destabilizing effect of the 
surface: two melting scenarios are found, surface melting 
(O2) with a continuous- and surface non-melting (Oi) 
with a discontinuous vanishing of the order parameter at 
the surface. More elaborate ab initio calculations reduce 
the problem to a mean-field order-parameter theorvand 
determine the free energy either as a lattice sum l7| or 
within a density functional theory (DFT) exploiting 
liquid-state correlations. 
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FIG. 1: Comparison of the melting line -Bm(r) obtained via 
the DFT-substrate approach (present work, full line) with 
that of Ref. [13 (dashed line). Top right: Numerical solutions 
of the order-parameter profile for different fields B increas- 
ing from top to bottom: (a) at T = 0.08 eorf (O2 transition 
with the liquid-solid interface invading the bulk as B y Bm) 
and (b) at T = 0.33 eod (Oi transition with the thick line 
corresponding to _B = Bm)- 



The vortex matter in Bi- and Tl- based high- Tc su- 
perconductors is characterized by an extreme anisotropy: 
pancake-vortices confined to superconducting layers ex- 
hibit repulsive logarithmic interactions within the planes 
and weakly attractive inter-layer interactions extending 
over many layers (see, e.g. Ref. 0; we ignore here a 
weak Josephson coupling between the layers). These 
anisotropic properties inspire the use of a substrate- 
mcan-ficld theory [lfij | describing the vortex system in 
terms of two-dimensional lattices of pancake- vortices sub- 
ject to a substrate potential generated by the mean ac- 
tion of the vortices residing in other layers. The bulk 
melting line is known 

[idi lUl llj to interpolate be- 
tween the Berezinskii-Kosterlitz-Thouless (BKT) transi- 
tion temperature Tqkt ^ of the individual layers at 
zero external field and the two-dimensional vortex-lattice 
melting temperature ^ Tc at high magnetic fields, 
see Fig. ^ We use the classical DFT of freezing |3| 
and exploit the anisotropic properties of the system with 
its natural separation of liquid-state correlations into in- 
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plane and out-of-plane components. As a result, we ob- 
tain a reliable and analytically tractable order-parameter 
theory which allows us to study (a6)-surface melting. We 
find two regimes within the B-T phase diagram: for low 
and high magnetic fields B the surface order parame- 
ter undergoes a finite, although reduced, jump (Oi) at 
the bulk melting temperature, whereas for intermediate 
values of B, stray magnetic fields destabilize the layers 
close to the surface, leading to the surface melting sce- 
nario (02)- 

We sketch the derivation of pancake-vortex interac- 
tions in a semi-infinite superconductor filling the half- 
space z > 0. These interactions are mediated by currents 
set up by the vortices via the Lorentz force (here, d is the 
layer separation, A the London penetration depth, and $0 
the flux unit): each vortex generates a current density 
j = — (c/47rA^)[(<l>o/27r)V</3 -I- A] which acts on another 
vortex core with a transverse force F = d^o z x j/c. For 
a co-planar pair of pancake-vortices the current density 
is driven by the phase gradient Wip — — x R/i?^ and 
the force oc 1 /R produces a long-range logarithmic repul- 
sion Vz^z{R) ~ — 2eo(iln(_R/^); this potential corresponds 
to that of a one-component charged plasma (OCP) with 
charge 2eod and £0 = ($o/47rA)^ the vortex hue 

energy. The interaction between two pancake- vortices re- 
siding in different layers derives from the vector potential 
A; calculating the field associated with a pancake- vortex 
in a semi-infinite geometry and integrating the Lorentz 
force provides the potential 
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with fz{K) = cxp(— _fir_^ 



and 

(3{K) = ~ K)/{K+ + K). The bulk term cx fz-z' is 
augmented by a stray-field term oc fz+z' relevant within a 
distance A from the surface. For small in-plane distances 
R ^ \, the contribution of the stray-field term can be 
neglected and we recover the bulk expression Vz^z'{R) w 



eod{dR/X^)[R/{R + 4\z- 



For a large separation 



i? 3> A, the surface term is relevant and we obtain 

Vz,z'iR) « 2eod[M^,z')ln{R/X) + (d/i?)e-(^+^')/^] , 

where 0t(^,z') = (d/2A)(e^l^-^'l/-^ + e-(^+^')A) < 1 is 
the fraction of flux trapped in the layer z' generated by 
a pancake- vortex at z: a test vortex at z' then effectively 
experiences a logarithmic attraction from two bulk-type 
pancake-vortices, the real one at z and a fake mirror vor- 
tex with equal sign at — z, the latter generated by the 
stray field. The algebraic repulsion associated with the 
second term in Vz.z' (R) is again due to the stray field and 
produces a surface softening. 

In our investigation of the vortex solid-liquid transition 
we make use of the classical density functional theory 



(DFT) [TjI which builds on the (grand canonical) free- 
energy difference SQ — figoi — ^^Uq expressed through the 
variation 5p{r) = p{r) — p in particle density away from 
the uniform liquid state density p; for the inhomogeneous 
and anisotropic vortex matter system, Spz(R) = p2(R) — 
p and p is the 2D liquid density, 
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|R-R'|)<5p,.(R') . (2) 



The first two terms describe the entropic contribution of 
a non-interacting gas, while the last term accounts for 
the microscopic interactions Vz^z' (R) via the direct pair- 
correlation function Cz,z'{R) of the liquid state; in the 
homogeneous liquid, the (dimensionless) Fourier trans- 
form Ck^iK) = (p/d) J d^r Cz{R)e^'^''^ is related to the 
static structure factor via <S'(k) = l/[f — Ck^{K)]. The 
appearance of finite density-modulations |(5pz(R.)| > 
in the solid is a consequence of these correlations. We 
exploit the anisotropy of the pancake-vortex system and 
separate Cz,z'iK) into in- and out-of-plane parts, 

Cz,z'{K) = c^°(if)d<5(z-z') - Vz,z'{K)/T, (3) 

where (K) denotes the correlation function of the 2D- 
OCP as obtained from standard Monte-Carlo simula- 
tions 01 . The second term accounts for the weak inter- 
plane interactions; within lowest-order perturbation the- 
ory [isj it is given by the Fourier transform of the out-of- 
plane interaction (jj), Vz,z'{K)/T = -a{K)[fz-z'[K) + 
(3{K)fz+z'{K)] with a{K) = 2neod^p/TX^K'^K+. 

We begin with the homogeneous bulk system. Con- 
sider the Fourier transforms 5pz(Kn)/ P for the density 
field and ^2(Kn) for the molecular field 'Cz(R) = 
ln[pz(R) / p], where Kn denote the reciprocal triangular- 
lattice vectors. For a correlation function c^°{K) decay- 
ing rapidly beyond the first reciprocal lattice vector, only 
the first two components need be retained and we can re- 
strict the Ansatz to the form SpzCR)/ p ~ riz+Pzg(R) and 
^^(R-) ~ Kz+S.zgi'R) withg(R) = X]|K„|=G<^^P(*^n-R-), 
G = i-K/^/Sa^ (oa the lattice constant, = 2$n/\/3_B). 
Furthermore, Kz and ^z are related to /x^ via 

Kz^-H^z), Pz=^'{^z)/6, (4) 

with <i>(^) — Ins"^/^ d'^Rexp[^g{R)] and $' its derivative 
(here, s denotes the 2D unit cell area). Inserting this 
Ansatz into the free-energy density lI2Jl, we obtain the 
reduced functional for a bulk system in the form (we use 
the normalization 5lu ~ Sfl/TSp with S the sample area) 
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with a = a(G) and fz-z' — fz-z'{G). The first term 
describes the free-energy density of individual layers 
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with in-plane correlations = c^°{G) and out-of- 
plane correlations described by the substrate potential 
-ajidz/d)f, = ~2a/dG+ with G+ ^ VCPTX^. 
Note that both and have to be understood as 
functions of the order-parameter fj.^ via (0J). The sec- 
ond non-local term in 6uj accounts for inhomogeneities of 
the order-parameter field /i^ and involves the dispersive 
'elastic coefficient' afz-z'- Here, we neglect the small 
change in density across the transition described by 77^; 
its inclusion involves a more elaborate analysis account- 
ing for the discrete nature of the particles [lj| . 

The bulk melting line B^{T) is obtained from minimiz- 
ing the functional 5uj [n] for a homogeneous order param- 
eter fjL. At large values of the inter-planar interaction 
is negligible and the melting temperature is given by the 
solid-liquid transition of the 2D-0CP as described 
by the free-energy density Q without substrate poten- 
tial (termed i5a;^°). For small ?° (large temperatures), 
6u^°{li) exhibits only the liquid minimum at fj-uq — 0, 
cf. Fig. 12 Lowering the temperature, the correlator ?° 
increases and SijJ^°{^) develops a second minimum at a 
finite value of fi describing the solid phase. At the critical 
value ?° = Cc = 0.856 the 'solid' minimum drops below 
the 'liquid' one and the high- field liquid-solid transition 
takes place. The comparison with numerical simulations 
[itI allows us to check the accuracy of our approach: 
Monte Carlo simulations [l3| show that the 2D-0CP 
freezes at « e^d/lO where the correlator assumes 
the value ?° ~ 0.77 < Cc] more sophisticated versions of 
DFT cure this discrepancy [l8| . 
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FIG. 2: Top left: direct pair-correlation functions at T /eod = 
0.08 for the 2D-0CP, c^"{K) from MC simulations and 
c^^h{K) = ^°{K) + 2a{K)/dK+ for the full 3D system at 
melting {B = _Bm). Bottom right: Free energy Suj^° {fi) /T 
for values = 0.80, 0.83, 0.845, 0.856 (= Cc critical value, 
thick line), 0.87. At melting the order parameter jumps from 

flsol ~ 0.51 to ^liq = 0. 

At low B {< B\ = $o/A^), the inter-planar cor- 
relations become relevant and the 2D correlations ?° 
are augmented by the substrate potential. The con- 
dition Cc = ?° + 2a/dG+ = ?° + (\/3eorf/27rT)[l -|- 



inyV3)B/Bx]- 

Bm{T) 

Bx 



^ yields the bulk melting line 
^/3eod 



V3 



27rT\cr 



- 1 



(6) 



Given the temperature dependence of the OCP correla- 
tor ?°(T) as an input, we obtain the full melting line as 
shown in Fig. ^ our result agrees well with those from 
numerical simulations |l2| and improves upon results ob- 
tained from other liquid-state closure schemes 11]. 

Surface melting involves a non-uniform order-para- 
meter field fiz defined in a half-infinite sample z > 0. 
The translation-invariant correlator oc fz-z' in now 
has to be replaced by the full expression Vz^z' {G)/T with 
additional mirror and surface terms. Minimization with 
respect to ^.z provides us with the integral equation 



^^.J^llh{^^z) 



r°° dz' - 

6a/ -^[fz^z- + fz+z']{iJ'z- ^^z') 



12aG 



G + G+Jo d 



'dz' 



fz+z' f^z' = 0. (7) 



The order parameter profile ^z can be calculated nu- 
merically and we show two typical examples in Fig. ^ 
while we find the surface order-parameter iiq to vanish 
smoothly at T = 0.08 eo^^ {O2 scenario), a residual jump 
in fiQ is obtained at the higher temperature T = 0.33 Sod 
(Oi). In the following, we locate the multi-critical point 
Tmc on the melting line B^^iT) separating the continuous 
O2- from the discontinuous Oi regime. 

To make progress analytically, we have to simplify the 
non-local terms in {Tj). Concentrating on the bulk, a 
gradient expansion of the term cx fz-z' leads us to the 
differential equation (we define dz^iz = Mz) 



ff,'^=d,jL0l^^i^^z)/T 



(8) 



with ^2 = (6a/d) / dz fzZ^ = 12a/ dGl- This local ap- 
proximation is valid if the profile fiz varies slowly over 
the extension G~^ of the kernel; we have verified that 
this condition is fulfilled for T > 0.04 epd (correspond- 
ing to i? < 0.5 B\). Equation ^ has to be com- 
pleted with a boundary condition. Close to Tmc, the 
order parameter fj,z is small near the surface; linearizing 
df^^5Lu'^°^/T « 6(1 — Cc)/iz in Q we obtain the equation 



2(1 + r) Jo 



dz' [fz 



z+z'l l^z 



(9) 



with r = dG+{l - cc)/2a and f3 = {G+ - G)/{G+ + G). 
While the solution of this type of integral equation is 
a non-trivial task in general, a straightforward solu- 
tion is possible in the present case due to the particu- 
lar exponential structure of the kernel. Taking the sec- 
ond derivative of ©, we obtain the differential equation 
= 6(1 — Cc)/X2, which shows no trace of the bound- 
ary term cx /? and thus coincides with the bulk equation 
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^ at small This equation then admits the ex- 

ponential solution A exp{^/rG^z) + B exp(— y^G+z) and 
inserting this Ansatz back into Q the boundary term 
fixes the ratio A/ B; as a result, we obtain the boundary 
condition 



K/Ai.].=o = G+(l-/3)/(l + /3)=G. 



(10) 



The analysis of the boundary value problem ||SJ) with 
(|10|l follows the one in Ref. 0|: Combining the bound- 
ary condition H10() with the expression for the 'conserved 
energy' {£ fi'J'^ - 2(5a;™^(^^)/T = 0, we find the relation 



(11) 



The liquid surface fiQ — is always a solution and we deal 
with a continuous surface melting (O2 scenario) if it is the 
only one. Once a second solution with /ip > is present, 
the surface undergoes a discontinuous Oi transition. The 
Oi- and O2 scenarios are separated by a multi-critical 
point: expanding Slo^^^, we find that admits two 
solutions for T > T^c ~ 0.29 eod (Bmc ~ 0.007Ba) and 
the surface (non)-melting is realized for T < Tmc {T > 
Tmc)- In Fig-El we show the solution of Hll|) and compare 
it with numerical results. 
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FIG. 3: Order parameter ^0 at the surface: analytic (see Eq. 
mi l, solid line) and full numerical solution (stars). Inset: en- 
ergy as a function of soliton position Zs- Note the minimum 
for T = 0.33 £orf > Tmc (Oi, pinned soliton) which has disap- 
peared at T = 0.28 eorf < Tmc (O2, depinned soliton). 

The appearance of the multi-critical point (T,„c, -B,nc) 
can be interpreted as a surface-depinning transition of the 
solid- liquid interface (soliton): the energy cost to push 
this soliton into the system generates a repulsive poten- 
tial, while the surface term favoring the liquid phase helps 
the soliton to enter, see Fig.|21 With decreasing temper- 
ature the stable minimum close to the surface (gener- 
ating a finite ^0 and leading to a Oi transition) moves 
deeper into the bulk and disappears altogether at Tmc- 
As a consequence, the (half-height) position Zs of the soli- 
tonic solid-liquid interface diverges logarithmically into 
the bulk as T ^ from below, Zs - | lii{Tm - T)\, 
within the entire O2 regime B > Smc. 

Finally, at high magnetic fields the layers melt indepen- 
dently following a first-order type 2D melting scenario. 



The order parameter /ig in the topmost layer then un- 
dergoes a finite jump and the surface non- melting (Oi) 
scenario applies, implying the existence of a second multi- 
critical point at high fields. Indeed, we do find such a 
jump at fields of order 10i?A, however, a more elaborate 
version of DFT is required for an accurate determi- 
nation of this second multi-critical point. 

In conclusion, we have analyzed surface melting in the 
pancake-vortex system of layered superconductors and 
have found both surface-melting (O2) and surface-non- 
melting (Oi) scenarios. The O2 scenario is realized over 
most of the low-field phase diagram and explains the ex- 
perimental observation of asymmetric hysteresis 0. 
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